##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--gene"), type = "character") ,
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--mut_rate_point_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    gene <- "RHOA"
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.Type.tsv",sep="")
    mut_rate_point_file <- paste(work_dir,"/images/mutRate/MutRate.RecurrentPoint.Type.tsv",sep="")
    info_file <- paste(work_dir,"/config/tumor_normal.class.list",sep="")
    images_path <- paste(work_dir,"/images/mutRatePlot/RHOA",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene <- opt$gene
info_file <- opt$info_file
mut_rate_gene_file <- opt$mut_rate_gene_file
mut_rate_point_file <- opt$mut_rate_point_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")

###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
dat_mutRatePoint <- data.frame(fread( mut_rate_point_file ))
dat_info <- data.frame(fread( info_file ))

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol==gene)
dat_mutRatePoint <- subset(dat_mutRatePoint , Hugo_Symbol==gene)
dat_mutRatePoint$Hugo_Symbol <- dat_mutRatePoint$vid
dat_mutRatePoint <- dat_mutRatePoint[,-3]

dat_plot <- rbind( dat_mutRateGene , dat_mutRatePoint )
dat_plot$Class <- factor( dat_plot$Class , levels = c("IM" , "IGC" , "DGC") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$MutRate , 2) * 100 , "%") 

dat_sampleNum <- dat_info %>%
group_by( Type ) %>%
summarize( SampleNum = length(unique(ID)) )
dat_sampleNum <- data.frame( dat_sampleNum )

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################
## 计算P值
if(gene=="AK2" | gene == "RHOA"){
    dat_plot <- dat_plot
}else{
    dat_plot <- subset( dat_plot , Hugo_Symbol %in% unique(dat_plot[dat_plot$MutNum >= 2,"Hugo_Symbol"]))
}

class_compare <- c("IM" , "IGC" , "DGC")

result <- c()
for(geneN in unique(dat_plot$Hugo_Symbol)){

    print(geneN)

    for( type in unique(dat_plot$Type)){

        dat_plot_tmp <- subset( dat_plot , Type == type & Hugo_Symbol == geneN )
        if(type == "IM + IGC"){
            class_compare <- c( "IM" , "IGC" )
        }else if(type == "IM + DGC"){
            class_compare <- c( "IM" , "DGC" )
        }else if(type == "IM + IGC + DGC"){
            class_compare <- c( "IM" , "IGC"  , "DGC")
        }
       # class_compare <- as.character(dat_mutRateGene$Class)

        if(nrow(dat_plot_tmp) > 0){
            for( i in 1:(length(class_compare)-1) ){
                class1 <- class_compare[i]
                for( j in (i+1):length(class_compare) ){
                    class2 <- class_compare[j]
                    tmp_1 <- subset( dat_plot_tmp , Class %in% class1 )
                    if(nrow(tmp_1)==0){
                        tmp_1 <- data.frame( Type = type , Hugo_Symbol = geneN , Class = class1 , MutNum = 0 , 
                            SampleNum = dat_sampleNum[dat_sampleNum$Type==type,"SampleNum"] , 
                            MutRate = 0 , value_text = "" )
                    }

                    tmp_2 <- subset( dat_plot_tmp , Class %in% class2 )
                    if(nrow(tmp_2)==0){
                        tmp_2 <- data.frame( Type = type , Hugo_Symbol = geneN , Class = class2 , MutNum = 0 , 
                            SampleNum = dat_sampleNum[dat_sampleNum$Type==type,"SampleNum"] , 
                            MutRate = 0 , value_text = "" )
                    }

                    tmp <- rbind( tmp_1 , tmp_2 )
                    tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
                    p <- fisher.test(tmp_fisher)$p.value
                    if( p < 0.001 ){
                        p_text <- trans(p)
                    }else{
                        p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
                    }

                    tmp$p.value <- p
                    tmp$p_text <- p_text
                    tmp$Type <- type
                    tmp$Hugo_Symbol <- geneN

                    result <- rbind( result , tmp )
                }
            }
        }
    }
}

result$Hugo_Symbol <- factor( result$Hugo_Symbol , 
    levels = unique(result[order(result$MutRate , decreasing=T),"Hugo_Symbol"]) , order = T)
result$Class <- factor( result$Class , levels = c("IM" , "IGC" , "DGC") , order = T )
out_name <- paste0( images_path , "/MutRate_" , gene , ".divide.tsv" )
write.table( result , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################

dat_use <- unique( result)
dat_use$Type <- factor( dat_use$Type , levels = c("IM + IGC" , "IM + DGC" , "IM + IGC + DGC") , order = T )
dat_use$p_text <- ifelse( dat_use$Type == "IM + IGC + DGC" , "", dat_use$p_text)
dat_use$p.value <- ifelse( dat_use$Type == "IM + IGC + DGC" , "", dat_use$p.value)
dat_use <- unique( dat_use)
dat_use <- subset( dat_use , Type != "IM + IGC + DGC" )

for( geneN in unique(dat_plot$Hugo_Symbol) ){

    dat_use_tmp <- subset( dat_use , Hugo_Symbol == geneN )
    dat_use_tmp$Type_num <- paste0( dat_use_tmp$Type , " (" , dat_use_tmp$SampleNum , ")" )
    plot <- ggplot( data = dat_use_tmp , aes( x = Class , y = MutRate , fill = Class ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Mutation Rate")+
    facet_grid(.~Type_num,space='free_x',scales='free_x') +
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1)+
    geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='none',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.ticks.x = element_blank(),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

    if(length(unique(dat_use_tmp$Type)) == 1 ){
        width = 2
    }else{
        width = 4
    }
    out_name <- paste0( images_path , "/MutRate_" , geneN , ".divide.pdf" )
    ggsave( out_name , plot , width = width , height = 5 )
}